#include "ppdef.h"

!--------------------------------------------------------
!----- subroutine readMatrix ----------------------------
!--------------------------------------------------------

! reads the data and mask files and creates 'matrix', 
! containing data and missing data, no land values
! 
!--------------------------------------------------------
subroutine ReadMatrix(sstfile,matrix,m,n,valex,maskfile,norma,imax,jmax,first,fileMean,fileStd)
  use ufileformat
  implicit none

  character(len=200),intent(in) :: sstfile(:),maskfile(:)
  real,pointer                  :: matrix(:,:)
  integer,intent(in)            :: norma
  integer,intent(out)           :: m,n
  real,intent(out)              :: valex
  integer, pointer              :: imax(:),jmax(:),first(:)
  real,pointer                  :: fileMean(:),fileStd(:)

  integer, allocatable          :: land(:),sea(:),nlines(:)
  real,pointer                  :: file(:,:,:),mask(:,:)
  real, parameter               :: valexc = 9999.
  integer                       :: i,j,s,t,k,d,q,nbvars,testN
  real                          :: Y,maskvalex,testmean
  real,pointer                  :: fileNorm(:,:,:),valexN



  interface
     subroutine stat(file,mask,fileNorm,mean,stddev)
       real,pointer               :: file(:,:,:),mask(:,:)
       real,pointer               :: fileNorm(:,:,:)
       real,intent(out)           :: mean,stddev
     end subroutine stat
  end interface


  nbvars=size(sstfile)

  allocate(imax(nbvars),jmax(nbvars),land(nbvars),sea(nbvars),nlines(nbvars),first(nbvars),fileMean(nbvars),fileStd(nbvars))


  do q=1,nbvars


     call uload(sstfile(q),file, valex)

     !test

     call uload(maskfile(q),mask,maskvalex)
     where (mask.eq.0) mask = valexc

     !where (spread(mask,3,size(file,3)).eq.valexc) file = valex

     do t=1,size(file,3)
        do j=1,size(file,2)
           do i=1,size(file,1)
              if (mask(i,j).eq.valexc) file(i,j,t) = valexc
           enddo
        enddo
     enddo

     !testN = count(file.ne.valex)
     !testmean = sum(file,file.ne.valex)/testN
     !write(stdout,*)'testMean',testmean,testN, first(q),nlines(q),maskfile(q),sstfile(q)

     !write(stdout,*)'testMean',sum(file,file.ne.valex

     !    ------------------------------------
     !     check if all input matrices
     !     have the same size as their masks
     !    ------------------------------------

     if(size(file,1).ne.size(mask,1).or.size(file,2).ne.size(mask,2)) then
        write(stdout,*) 'File', trim(sstfile(q)), 'does not have the same size as its mask,', trim(maskfile(q)), 'please check.'
        STOP
     end if

     imax(q)=size(mask,1)
     jmax(q)=size(mask,2)

     !    -----------------------------------------
     !     temporal dimension, same for all files!!
     !    -----------------------------------------  
     nlines(q)=size(file,3)


     !    -----------------------------------------
     !     counts land points
     !    -----------------------------------------    
     land(q)=count(mask.eq.valexc)

     !    -----------------------------------------
     !     calculates first dimension of the matrix 
     !     (spatial points)
     !    -----------------------------------------    
     sea(q)=(imax(q)*jmax(q))-land(q)



     if(q.eq.1) then
        first(q)=1
     else
        first(q)=first(q-1)+sea(q-1)
     end if

     deallocate(file,mask)

  end do

  if(any(nlines.ne.nlines(1))) then
     write(stdout,*) 'The temporal sizes of both files provided are not the same'
     write(stdout,*) 'I have found the following temporal dimensions,', nlines(:)
     write(stdout,*) 'Please bear in mind that this routine is constructed for the'
     write(stdout,*) 'case when the temporal size of all files being reconstructed is the same.'
     write(stdout,*) 'If you wish to do more complicated things, please do so, but now, put in'
     write(stdout,*) 'two files with the same temporal size (and the same temporal distribution'
     write(stdout,*) 'if you want the results to be meaningful) Thanks.'

     STOP
  end if




  write(stdout,*)
  write(stdout,*) '****************************************************************'
  write(stdout,*) 'Now some statistics about your data:'
  write(stdout,*)
  write(stdout,*)
  write(stdout,*)'number of land points in the mask', land(:)
  write(stdout,*)
  write(stdout,*)'dimension of the files entered (1st spatial dimension)'
  write(stdout,*) imax(:)
  write(stdout,*)  
  write(stdout,*)'dimension of the files entered (2nd spatial dimension)'
  write(stdout,*) jmax(:)
  write(stdout,*)
  write(stdout,*)'dimension of the files entered (temporal dimension)   '
  write(stdout,*) nlines(:)
  write(stdout,*)

  m=sum(sea(:))
  allocate(matrix(m,nlines(1)))


  do q=1,nbvars


     call uload(sstfile(q),file, valex)
     where (file.eq.valex) file = valexc 


     call uload(maskfile(q),mask,maskvalex)
     where (mask.eq.0) mask = valexc

     !    -----------------------------------
     !     normalisation of the input matrices
     !     if norma = 1
     !    ------------------------------------

     if(norma.eq.1) then
        call stat(file,mask,fileNorm,fileMean(q),fileStd(q))
        !fileNorm=file
     else
        call stat(file,mask,fileNorm,fileMean(q),fileStd(q))
        fileNorm=file
     end if


     where(fileNorm.eq.9999) fileNorm=valexc

     ! where (spread(mask,3,nlines(q)).eq.valexc) fileNorm = valexc

     ! testN = count(fileNorm.ne.valexc)
     ! testmean = sum(fileNorm,fileNorm.ne.valexc)/testN

     !    ---------------------------------------
     !     reshape input matrices into one matrix
     !     size M*N (M spatial dimension,
     !               N temporal dimension)
     !    ---------------------------------------

     do t=1,nlines(q)

        s = first(q)

        do i=1,imax(q)
           do j=1,jmax(q)
              if (mask(i,j).eq.valexc) fileNorm(i,j,t) = valexc

              if(mask(i,j).eq.valexc.and.fileNorm(i,j,t).ne.valexc) then
                 write(stdout,*) 'WARNING ',i,j
              end if



              if(mask(i,j).ne.valexc) then

                 matrix(s,t)=fileNorm(i,j,t)

                 s = s+1

              endif
           enddo
        enddo
     enddo


     deallocate(file,fileNorm,mask)

  enddo

  valex=valexc

  n=nlines(1)


  !    -----------------------------------
  !     write mean and standard deviation
  !     of input matrices
  !    ------------------------------------

  write(stdout,*)
  write(stdout,*)'mean',fileMean(:)
  write(stdout,*)'Standard deviation',fileStd(:)
  write(stdout,*)


deallocate(land,sea,nlines)

end subroutine ReadMatrix





